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c/2\,BSTRACT: We used computer simulations to study spontaneous strain localization in granular materi- 

'^Is, as a result of symmetry breaking non- homogeneous deformations. Axisymmetric triaxial shear tests 

rtvere simulated by means of standard three-dimensional Distinct Element Method (DEM) with spheri- 

\taX grains. Carefully prepared dense specimens were compressed between two platens and, in order to 

imic the experimental conditions, stress controlled, (initially) axisymmetric boundary conditions were 

onstructed. Strain localization gave rise to visible shear bands, previously found experimentally under 

C^imilar conditions by several groups, and different morphologies could be reproduced. We examined 

' rtie stress-strain relation during the process and found good agreement with experiments. Formation 

,-^echanism of shear bands is discussed. 
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INTRODUCTION 

^Spontaneous symmetry breaking in granular mate- 
^>Oials occurs in many different forms. Here we focus 
^)n strain localization and subsequent development 
'^f shear bands. Shear bands appear nearly always 
rtf dry granular material is subjected to shear. Its 
iHrst study dates back to the nineteenth century 
'"^nd since then it was investigated in many differ- 
ent geometries and specially designed laboratory 
Cjiests (e.g. plane strain, biaxial, and triaxial tests). 
^jHere we present numerical studies of axisymmetric 
"^riaxial tests, the most common laboratory tests in 
''^_G(eomechanics. 

. 5^ I In a simplified picture, a triaxial test typically 
consists of a cylindrical specimen enclosed between 
two end platens and surrounded by a rubber mem- 
brane. An external pressure is applied on the mem- 
brane, either by placing the system into a pressur- 
ized fluid, or creating a relative vacuum inside the 
system. The end platens are pressed against each 
other in a controlled way, either with constant ve- 
locity (strain control) or with constant force (stress 
control). The force resulting on the platens, or the 
displacement rate of the platens is recorded, as well 



as the volume change of the specimen. 

The triaxial test is an elementary test, per- 
formed to obtain mechanical properties of soils. 
Antifriction devices (lubricated end platens) were 
designed in order to suppress strong heteroge- 
neous responses, such as barreling and local- 
ization of deformation along failure planes. In 
the past 20 years the study of localization pat- 
terns gained more attention and strain localiza- 
tion became an important research field, as ex- 
perimental tools as Computed Tomography (CT) 
became available to study the internal struc- 
ture of strained s pecimens (jDesrues et al. 1996| 
IBatiste et al. 2004|) . Such studies revealed complex 
localization patterns and shear band morphologies 
depending on the test conditions. 

2 SIMULATION METHOD 

We used standar d th ree-dimensio nal D istinct 
Element Method (fCun dall Strack T979|l (also 
known as Molecular Dynamics) to perform sim- 
ulations of strain controlled triaxial shear tests. 
As an advantage, contrary to the Finite Ele- 
ment Method (FEM) commonly used in sim- 
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Figure 1. Simulation setup. A cylindrical specimen placed 
between rigid horizontal platens and surrounded by an elas- 
tic membrane is subjected to axial load and confining pres- 
sure. 

ulations of soils, the Distinct Element Method 
(DEM) does not require any macroscopic constitu- 
tive model, instead it is based on a microscopic con- 
tact model (for a review see ( Luding 2004D and ref- 
erences therein). We u sed the Hertz contact model 
fILandau fc Lifshitz 19701 with appropriate damp- 
ing fBrillia ntov et al. 1996)) . combined with a fric- 
tional spring-dashpot model ( Luding 20 04). 

The simulation setup can be seen on Fig. ^ An 
initially axisymmetric system of spherical grains 
(particles) was placed between two horizontal 
platens. The bottom platen was fixed. On the up- 
per platen, having mass M = 10~^ kg, an axial 
load was applied (as described later). The upper 
platen could not rotate along the vertical axis. The 
rotational inertia of the upper platen in tilting was 
/ = 10~^ kg m? . The particles were surrounded 
by an "elastic membrane" composed of overlapping 
spheres having equal diameter dm = 10~^ m and 
mass density pm = 0.1 ■ 10^ kg/m?, and initially 
forming a triangular lattice (see part (b) of Fig.^). 
The rotational degree of freedom of the "membrane 
nodes" was frozen (i.e. they could not rotate). 

The "membrane nodes" were interconnected 
with linear springs having zero base length and ini- 
tial elongation /q = 0.5-10"^ m (equal to the initial 
distance of neighboring "membrane nodes"). The 
force Fs, acting between two "membrane nodes" 
connected with a spring, was calculated as 



(1) 



where Hg = 0.5 N/m and 7^ = 10^^ N s/m 
are stiffness and damping coefficients, Ig is the 
spring's elongation and Vs is the relative veloc- 
ity of the nodes. At any time the spring's elonga- 
tion is equal to the relative distance of the nodes. 
The stiffness of the springs was chosen such that 
the particles could not escape by passing through 
the membrane. Additionally a confining pressure 
(Tc = 0.5 ■ 10^ N/m? was applied on the membrane, 
by calculating the forces acting on the triangular 



facets formed by neighboring "membrane nodes 
(see part (c) of Fig. H}. The effective external pres- 
sure in this setup is larger than ctc, because the 
"membrane springs" have their own contribution 
as well, however, with the used parameters, this 
contribution is small compared to Uc. 

For simplicity we made no difference between 
particle-particle, particle-platen, and particle- 
membrane contacts. The normal F„ and tangential 
Ft components of the contact force were calculated 
as 



(2) 
(3) 



where k„ = 10^ N/m^/"^, Kt = 10^ N/m, 7„ = 
1 s/rrc'l'^ , and 74 = 1 N s/m are the normal 
and tangential stiffness and damping coefficients, 
5n and 6t are normal and tangential displacements, 
and Vn and are the normal and tangential rel- 
ative velocities. The normal displacement 5„, the 
normal velocity f„, and the normal force F^ are 
one-dimensional quantities measured along to the 
normal vector of the contact plane, while the tan- 
gential displacement 5t, the tangential velocity v^, 
and the tangential force are two-dimensional 
vectors embedded in the contact plane. 

Knowing the relative position, the shape, and 
the size of the bodies in contact, the normal dis- 
placement can be calculated directly. Calculating 
the tangential displacement is much more compli- 
cated: The tangential velocity must be integrated 
during the lifetime of the contact. This integration 
must be performed in the contact plane. In our 
simulation program the tangential displacement 5t 
was implemented clS Qb 3D vector in the observa- 
tional space. During integration, this vector was 
rotated as the local configuration changed, keep- 
ing it always in the contact plane. 

The Coulomb friction law limits the frictional 
force to pFn, where p = 0.5 is the coefficient 
of friction. To allow for sliding contacts, we also 
limited the length of the tangential displacement 
to fiFn/ Kt, shortening the displacement vector ac- 
cordingly. This frictional spring-dashpot model im- 
plements both sliding and static friction. 

Our simulations are run at zero gravity. After 
calculating the interaction forces, and adding the 
external load and confining pressure, the motion 
of bodies (grains, "membrane nodes", and upper 
platen) is calculated by solving numerically the 
Newton equations using a given At = 10~^ s in- 
tegration time step. The translational motion is 
calculated with Verlet's leap-frog method. The ro- 
tational state of bodies (given in quaternion rep- 
resentation) is integrated with Euler's method. 
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In quasi-static processes as the one simulated by 
us, the vibration introduced by grain (spring) elas- 
ticity is basically undesired noise. We checked the 
noise level in our simulations and set the param- 
eters to keep it low. The inverse of the eigenfre- 
quency of all contacts, in both normal and tangen- 
tial direction, is more than one order of magnitude 
larger than the integration time step. 

In the first part of the simulation (prepara- 
tion phase) a hard cylinder touching the internal 
side of the membrane was introduced, and thus 
the membrane was neglected. The particle system 
was built by randomly placing spheres into this 
cylinder. The maximum allowed initial grain over- 
lap was 1%. To assure reasonable execution time, 
we started with a sufficiently tall system (usu- 
ally 3 time taller than the final system size). The 
particles (grains) were given equal mass density 
pp = 7.5 ■ 10^ kg/m?. The diameter of the spheri- 
cal grains was taken from a Gaussian distribution 
with mean value (rf) = 0.9 ■ 10~^ m and standard 
deviation Ad = 0.025 ■ 10"^ m, and cut at 4A(i 
around the mean value. 

In the preparation phase the rotational degree 
of freedom of the upper platen was frozen and the 
coefficient of friction was set to zero. Each particle 
and the upper platen were given a velocity propor- 
tional to a contraction velocity f c = 80 ■ 10~^ m/s 
and their distance from the bottom. The system 
contracted until the upper platen's velocity in 
grain-platen collisions decreased to zero. At this 
point an axial load Fq = 200 ■ 10~^ was switched 
on, which further contracted the system. After the 
system relaxed, the cylinder was removed, letting 
the membrane and the confining pressure carry the 
load. The axial load Fq and the confining pressure 
(Tc were chosen such that the system came to equi- 
librium without barreling. We prepared one sam- 
ple having diameter D = 22 ■ 10~^ m and height 
H = 46 ■ 10~^ m, containing Np = 27000 particles 
and = 14904 membrane nodes. The sample's 
geometry factor H/D ^ 2, is similar to the typical 
geometry factors used in experiments. 

For preparing sphere packings 
(jWeitz 2004) many different methods ex- 
ist (e.g. dLubachevsky fc Stillinger 1990| 

ISherwood 1997|) ). Using the deposition method 
described above, the volume fraction at the end 
of the preparation phase was /o = 0.643, which 
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Figure 2. Stress-strain relation. The stress ratio (ct/cto, 
where ctq denotes the initial stress) measured on the up- 
per platen is shown as function of the axial strain, for the 
executed simulation runs (see Tab. (See inset for low 
stress ratios.) For the lower two curves (a, c) tilting of the 
upper plate was enabled, and for the upper two (b, d) it 
was disabled. 

is slightly larger than the random close packing 
value of identical spheres, as expected for an 
ensemble of spheres with size dispersion. 

After the preparation phase the sample was 
compressed by moving the upper platen down- 
ward in vertical direction with constant velocity 
Uc (strain controlled experiment). In the executed 
simulation runs we used two different velocities: 
Uc = uq = 10~^ m/s (base value), and Uc = 2uq 
(two times faster). During compression, tilting of 
the upper platen was either enabled or disabled. 
We executed four different runs (see Tab. H]) start- 
ing from the same initial condition. 

During the runs we measured the stress a on the 
upper platen, and calculated the stress ratio cr /aQ, 
where (Jq denotes the initial stress. The identifica- 
tion of the shear bands is a non-trivial task. One 
possibility i s to generalize the m ethod used in two 
dimensions (jDaudon et al. 1997| . i.e., to calculate 
the shear intensity around the particles from the 
local deformation tensor. After doing this, we have 
realized that monitoring the rotational state of the 
particles is sufficient for the purpose of shear band 
identificiation, and it leads to the same results. Our 
observations are presented in the next section. 

3 SIMULATION RESULTS 

As the axial strain increases, the response of the 
granular sample (the stress ratio) increases until it 
reaches a peak value (see Fig. Ej). This is a basic 
observation of triaxial shear tests of dense gran- 
ular specimens (Lade 2002). Strain localization in 
granular materials is followed by a decrease in load 
bearing capacity. Dense granular materials dilate 
during shear. Due to this, the load bearing parti- 
cle chains collapse. Shear bands occur after peak 
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Figure 3. Vertical cross sections. The snapshots were taken 
at the middle of the sample at 10% axial strain from the 
(c) (left) and the (d) (right) simulations. The darkness (red 
content of color) is proportional to the rotational energy. 

failure and result in further decrease in strength. 

According to Fig. 121 up to 15% axial strain there 
is no significant difference in the stress-strain re- 
lation measured in the different simulation runs, 
indifferent of the strain rate and tilting of the up- 
per platen. However, we observe a change after this 
point. At large axial strain values, when tilting of 
the upper platen was disabled, the stress ratio in- 
creases again. 

Taking cross sections of the sheared samples, 
and coloring the grains according to their rota- 
tional state we could visualize the shear bands (see 
Fig. 121). On these renderings the larger the rota- 
tional energy, the more red (or darker in black-and- 
white version) a particle gets. The fact that the 
rotational state of the grains identifies the shear 
bands, is due to the fact that the shear bands are 
characterized not only by dilation but also by ro- 
tation of the grains, which is known from both ex- 
periments and simulations (jOda fc Kazama 1998| 
Herrmann et al. 2nn"4|) . 

Another important result is that internal in- 
stabilities can develop into a symmetry-breaking 
localized deformation along a failure plane when 
tilting of the upper platen is enabled, while 
nontilting platens act stabilizing factor 

resulting in two axisymmetric conical surfaces 
and complex localization patterns around them 
fjFazekas et al. 2005|) . This result is confirmed by 
similar experiments executed in micro-gravity 
fIBatiste et al. 2nn4|) . and proves that in the ab- 
sence of reinforced axisymmetry, spontaneous sym- 
metry breaking can take place as a result of inter- 
nal instabilities (jPesrues et al. 1996|1 . 

4 CONCLUSIONS 

We executed triaxial shear test simulations using 
DEM. Different shear band morphologies known 
from experiments could be reproduced. To our 



knowledge it is the first time that these localization 
patterns were reproduced in DEM simulations. We 
showed that the shear bands can be identified with 
the rotational state of the grains, and symmetry 
breaking strain localization can develop if the sym- 
metry is not enforced with nontilting platens. The 
agreement of our results with the experimental re- 
sults is very good, even if the system size (num- 
ber of particles) in our simulations is much smaller 
than in experiments. 
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